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ABSTRACT 


A computerized solution for three dimensional stress 
analysis was developed. The solution was based on the finite 
element technique employing twenty and thirty-two nodal 
point three dimensional isoparametric elements. The solution 
is applicable to linear elastic structures composed of homo¬ 
geneous and isotropic materials subjected to static loading 
under constant temperature. Deformed states are restricted 
to small displacements and displacement gradients. The 
program is coded in standard FORTRAN IV and is machine inde¬ 
pendent except for one subroutine which may be adjusted to 
specific installations. „ Independence of core space require¬ 
ments is achieved by use of direct access storage facilities. 
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I. INTRODUCTION 


The finite element technique is a relatively new method 
for stress analysis of structural continua. 

Due to its almost limitless abilities this technique has 
recently been extensively used in structural analysis and 
considerable research has been conducted toward the develop¬ 
ment of elements and techniques satisfying the requirements 
of advanced structural technology. 

The present work concerns the development of a general 
purpose digital computer program for three dimensional stress 
analysis, named TRISOP, and employing the finite element 
method. 

The program developed, in its present form, is applicable 
to the solution of problems concerning linearly elastic struc¬ 
tures, composed of homogeneous and isotropic materials, and 
subjected to static loading at constant temperature. The 
deformed states are restricted to small displacements and 
displacement gradients. The use of an equation solver of, in 
principle, infinite capacity devised by Cantin [Ref. 2] in 
conjunction with direct access storage devices enabled the 
memory space requirements to be independent of the size of 
the problem, the latter being limited only by the execution 
time and availability of storage devices. 

At the time of initiation of the development of TRISOP 
the most promising finite elements, in the three dimensional 
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analysis field, were the family of isoparametric elements. 

Two such elements, namely the 20 (Quadratic) and 32 (Cubic) 
nodal points elements, are used in the program. 

The present form of TRISOP is an improved version of a 
program coded by Professor Gilles Cantin of the Naval Post¬ 
graduate School and employing quadratic isoparametric finite 
elements. The major improvements brought to this initial 
form of TRISOP by the present version are, the extension to 
cubic elements, inclusion of predefined displacements in the 
boundary conditions, improved execution time and independence 
from the size of the problem. 

The structure of the program follows a modular pattern 
in order to provide for further improvements and, possibly, 
adjustment to new developments in the field. 
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II. NATURE AND FORMULATION OF THE PROBLEM 


The problem considered is a three dimensional stress anal 
ysis of a statically loaded structure under constant tempera¬ 
ture. The analysis is restricted to deformed states within 
the elastic region and involving small displacement gradients 
Consequently the linear strain tensor only need be considered 
and the strain-displacement relations can be expressed as: 
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( 2 . 1 ) 


or in more compact form: 

{e} = [D]{u} 


( 2 . 2 ) 


Where both {c}iand {u} are point functions. 

As a consequence of the above the principle of super¬ 
position applies to strains and displacements. The stress- 
strain law is assumed to be linear and thus the principle of 
superposition applies for stresses as well. Therefore, 
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displacements, strains and stresses of two distinct deformed 
states are algebraically additive resulting into a new 
compatible deformed state. Consequently initial stresses and 
strains need not be considered in the formulation of the 
problem. If such strains and stresses exist they will have 
to be superposed to the results of the solution which will be 
relative to a reference state, namely the state of the struc¬ 
ture before the application of the load considered in the 
problem. 

In addition to small displacement gradients, the analysis 
is restricted to displacements small enough that equilibrium 
and compatibility relations can be referred to the undeformed 
geometry of the structure (here and in what follows by unde¬ 
formed state is meant the reference state) . 

Due to the restriction of small displacements the question 
of stability does not enter the analysis [Ref. 7]. The struc¬ 
ture may be composed of more than one materials which are 
assumed to be homogeneous and isotropic. The constitutive 
laws for this problem can be expressed in matrix form as: 
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or in more compact form; 


{a} = [E]{e} (2.3) 

Where {a} is a point function, E* is Young's modulus and \i 
is Poisson's ratio. 

A. FINITE ELEMENTS 

The main feature of the finite element technique is the 
approximation of the elastic continuum possessing an infinite 
number of degrees of freedom by an assemblage of substructures, 
the finite elements, each possessing a finite number of degrees 
of freedom. 

Among the various types, two elements of the isoparametric 
family were chosen for this development; namely, the 20 and 
32 nodal point three dimensional elements. The concepts under¬ 
lying the construction and use of these elements as well as 
the resulting properties will not be presented here except for 
those necessary for the presentation of the structure of the 
computer program. An extensive treatment of this family of 
finite elements as v/ell as the isoparametric concept may be 
found in Refs. 4 and 8. 

Two systems of reference are used throughout. First, a 
"local" non dimensional and, in general, non cartesian system 
is used for convenience of computation of the element 
stiffness matrices. Then everything is transferred in a 
"global" cartesian system of reference (x,y,z) for the final 
computations of the assembled system. 
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A numbering convention for the element nodes is employed 
and is shown in Figures 1 and 2. A different convention is 
used for numbering the nodes of the assembled structure and 
is described in Section IV.A. 

By use of shape functions describing the varia¬ 

tion of element properties in the local system, every element 
property 4>(5/n/C) can be expressed in terms of its nodal 
values as: 


Thus the displacements of points within an element will be 
given by: 
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or in more compact form: 

{u} = tN]{u^} 

An extensive treatment of the shape function concept may be 
found in Ref. 8. 
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The shape functions used in this development belong to the 
Serendipity family and are shown in Refs. 4 and 8. 

B. ELEIMENT AND TOTAL STIFFNESS MATRICES - CONSISTENT LOADS 
Every individual finite element can be visualized as an 
independent structure and thus the derivation of a stiffness 
matrix is possible. 

The derivation of the element stiffness matrix [K 1 is 

e 

shown in Ref. 1 where the following integral expression 
is finally obtained: 

1 1 1 T 

[K ] = / / / [N*]^[E][N*] det[J]d^dndc (2.4) 

-1 -1 -1 

where: [N*] = [D] [N] 

and [J] is the Jacobian matrix of the transformation. 

c 

This expression for the element stiffness matrix is ideally 
suited for Gaussian numerical integration. Computational con¬ 
siderations on the performance of this integration as well as 
the choice of integration points are presented in Ref. 1. 

After formation of the stiffness matrix for each element, 
the total stiffness matrix [K] for the structure can be formed 
by superposition of corresponding stiffness coefficients. 

The force-displacement relations for the whole structure 
can then be expressed as: 

{F} = [K]{u} (2.5) 

where{u}is a vector containing all nodal displacements of the 
structure and {F} >a vector of consistent concentrated loads 
applied at the nodes and equivalent to the loading applied to 
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the structure. The determination of the consistent loads is 
not incorporated in the present solution and it is assumed 
that it is known. 

Relations (2.5) represent a system of linear equations. 

When this system is initially formed, both {F} and {u} will in 
general contain known and unknown elements. By further 
manipulation the system can be brought to the standard form 
with all unknowns on one side and its solution will result 
into the determination of all unknown nodal displacements and 
forces (see Section III.B.5). 

C. STRESSES AND STRAINS 

By use of the strain-displacement relations (2.2) and the 
displacements obtained from the solution of system (2.5) the 
vectors of nodal strains can be determined. 

However the application of relations (2.2) is only possible 
at the element level in the finite element theory. Thus the 
nodal strains for each element separately are first evaluated. 

Further by applying the stress-strain relations (2.3) the 
element nodal strains are evaluated. 

Except for nodes lying on the interfaces between elements 
of different materials, a stress analysis requires the deter¬ 
mination of average nodal strains and stresses of all elements 
sharing common nodes. This is accomplished by simply averaging 
the algebraic values of the stress and strain components at 
each node. 
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III. DESCRIPTION OF THE COMPUTER PROGRAM 


Computerized solutions to structural problems based on 
the finite element technique have, in the recent years, led 
to a considerable amount of research in three disciplines: 
structural analysis, numerical analysis and computer program¬ 
ing. Thus new elements and computational techniques are 
continuously being developed paralleled by the development of 
computing equipment. For this reason it is not unusual for 
such computer programs to become obsolete within a few years 
of their completion. 

Under these circumstances, the development of a successful 
program is required to follow rules which will preserve its 
effectiveness as long as possible by an adjustment process. 
These rules are sxammarized to: 

a. Standard coding and machine independence enabling 
the program to be used by various installations 
regardless of the computing equipment available, and 

b. Modular construction allowing adjustment to new 
developments without requiring the construction of 
the entire program. 

These principles have governed the development of TRISOP 
at all stages that brought it to its present form. 

The structure of TRISOP follows closely the analytic devel¬ 
opment of Section II. (Also see Flow Diagram.) 
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The program is coded in standard FORTRAN and is machine 
independent except for the subroutine DISK via which transfer 
of data to and from direct access storage devices is accom¬ 
plished. All calculations are performed in double precision 
in order to provide adequate accuracy in spite of the round 
off errors. 

The memory space requirements are independent of the size 
of the problem, the latter being limited only by the size of 
direct access storage space available and execution time. 

The extension of the capabilities of the program, as for 
example inclusion of thermal stress analysis and nonhomoge- 
neous materials will easily be accomplished without inter¬ 
fering with parts of the program that are not altered by these 
changes. 

A number of checking points was incorporated in order to 
assure the proper development of the execution. Whenever 
required, error and/or warning messages are given in the out¬ 
put. If the continuation of the solution is useless due to 
the error detected the execution for this problem is stopped 
and the solution of the next problem (if such exists) begins. 

The most important variables used throughout the program 
are shown in Appendix A and a brief description of their 
function is given. 

The basic characteristics of the two programs, for 20 and 
32 nodal point elements, are the same. The minor existing 
differences are discussed in Appendix G where listings of both 
programs are given. 
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In what follows the function of the different modules of 
the program is explained with reference to the analytic devel¬ 
opment of the solution as it is presented in Section II. 

A. CONTROL PROGRAM 

The control program serves no other purpose but to initiate 
the calling sequence of the different modules of the program 
and read, from input data cards, the characteristics (number 
and length of records) of the data files as these are 
requested by the JCL^ cards. These characteristics are later 
on (in subroutine INPUT) checked for correctness and adequacy. 

B. SUBROUTINE MAINE 

Subroutine MAINE is the actual main program. It is coded 
in subroutine form in order to provide for the solution of 
more than one problem in one run by simply restarting the 
calling sequence of the different modules. Its function con¬ 
sists of a step by step calling procedure. The only numerical 
operation executed by subroutine MAINE itself is the computa¬ 
tion of the average nodal stresses and strains after the element 
nodal stresses and strains have been computed by subroutine 
STRESS. In the coding of subroutine MAINE a timing procedure 
has been incorporated giving at the output the execution time 
for each step. This is an optional feature not connected with 
the execution of the problem itself and can be removed if re¬ 
quired. However, it provides useful data for studying the 
factors affecting the execution time of a problem. 


JCL = Job Control Language 
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The calling procedure and the function of the subroutines 
called are as follows: 

1. Reading of Input Data 

/ 

By calling subroutine INPUT all input data cards are 
read. First, the general characteristics of the problem (as 
number of elements, number of nodal points, etc.) are read 
and used for the computation of other characteristics (as 
total number of equations, number of blocks per column). All 
general characteristics of the problem are stored in a common 
storage area and used throughout the program by most of the 
subroutines. Subsequently the arrays of connectivity, nodal 
point coordinates, material properties, concentrated load data, 
and boundary conditions are read. All arrays except for the 
material properties (which is of pre-specified size) are stored 
in data files on direct access storage devices. This feature 
enables the program to be independent of the size of the 
problem. The bandwidth and the number of blocks per row of 
the total stiffness matrix are computed from connectivity data. 
All input data as well as those computed are printed out for 
checking purposes. 

The nature of the program allows the solution of a great 
variety of problems. For this reason and in order to avoid 
limitations on the form of the input data most reading formats 
are not specified but left to the choice of the user by use 
of object time formats. 

Within subroutine INPUT several checks are made in order 
to insure the consistency of some data. All such checks 
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throughout the program are accomplished by calling the proper 
entry of subroutine ERROR. The checks made by subrouting 
INPUT are the following. 

a. Stop Trap for the Last Problem 

The last problem of each run is followed by an 
indicator card, as explained in Section IV. Once this indi¬ 
cator card is read, the executi.cn is stopped. 

b. Number of Different Materials 

The maximum permissible number of different mate¬ 
rials present in the structure to be analyzed is ten. If a 
problem involves a number of materials greater than ten, the 
execution is stopped and an error message is given in the 
output. 

c. Characteristics of Data Files 

The characteristics of data files, namely the 
number and size of records for each file, are given to the 
program once for every run in the main program but their con¬ 
sistency is checked for every problem of the run. If the size 
of one or more of these characteristics is smaller than that 
required for the problem the execution is stopped for this 
problem and execution of the next problem (if such exists) 
begins. In case the size of one of the data files character¬ 
istics exceeds that required for the problem, a warning message 
is given at the output and execution continues. 

2. Formation of Element Nodal Coordinate Arrays 

The coordinates of nodal points are provided to the 
program via subroutine INPUT in a sequence following the node 
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numbering convention for the whole structure. However, at 
various points of the program the use of arrays containing 
the nodal coordinates of each individual element is required. 
The formation of such arrays is accomplished by subroutine 
SORT. The sequence of elements of these arrays follows the 
numbering convention of the element nodes. Each array formed 
is stored in a direct access data file. 

3. Formation of Element Stiffness Matrices 

The individual element stiffness matrices [Ke] are 
formed by subroutine FSTF in the following way. 

The element stiffness matrix is given by the 
expression 

111 T, 

[Ke] = / / / [N*] [E] [N*] det [JjdCdndi; (3.1) 

-1 -1 -1 

First the matrix of elastic constants [E] is computed 
by calling subroutine ELPROP. In case the whole structure is 
made of one material this operation is executed only once 
since all elements have a common elastic matrix. In the case 
of more than one material matrix [E] is computed every time 
the element considered is of different material than the 
previous one. Subsequently subroutine CUB is called to per¬ 
form Gaussian numerical integration. For each Gauss point 
subroutine FORMK is called by CUB and the integrant of 
expression (3.1) is formed in a step by step procedure and 
evaluated at the integration point. Once an element stiffness 
matrix is formed it is stored on disk. 
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4. Formation of Total Stiffness Matrix 


The formation of the total stiffness matrix is per¬ 
formed in subroutine MERGE by superposition of the correspond¬ 
ing stiffness coefficients. In order to avoid the presence of 
the whole stiffness matrix in core, it is partitioned in blocks 
of pre-specified size. These blocks are formed one at a time 
by selecting the proper elements from the element stiffness 
matrices. Each element stiffness matrix is scanned in order 
to detect elements belonging to the block under formation. 

Every completed block of the total stiffness matrix is in 
turn stored on disk. The method of partitioning and storing 
of the total stiffness matrix is established in Ref. 2. A 
thorough study of this paper is essential in order to compre¬ 
hend the procedure outlined here as well as to become familiar 
with the program itself. 

5. Formation of Load Vector and Application of Boundary 

Conditions 

After the formation of the total stiffness matrix the 
completion of system 2.5 requires the formation of vectors {F} 
and {u} and some further manipulation of the system in order 
to bring it in the standard form. 

The load vector {F} is partitioned into a number of 
smaller vectors the size of which corresponds to the size of 
total stiffness matrix blocks. Each such vector is formed in 
sequence by use of the consistent load data, stored on disk 
from the first step of the execution (subrouting INPUT), and 
subsequently is stored on disk. Load vector elements 
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corresponding to unknown forces, namely reactions at the 
constraints, are replaced at this stage by zeros. 

The boundary conditions imposed to the problem, result 
in the presence of a number of unknown reactions, designated 
hereby F^*, and predefined displacements, designated by 
If such a boundary condition is applied to the ith element of 
the force and displacement vectors the system will have the 
fom; 
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This system is equivalent to a good degree of approximation 
to the system 
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Where B is a number of several orders of magnitude greater 

Fi* 

than the greater stiffness coefficient and F^** = - ^ 

The products ,...,F^**K^^ are negligible in size com¬ 

pared to the other coefficients. After the system is solved 
the unknown reaction F^* is determined by: F^* = - BF^** 
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This is a well known method of solution of mixed systems 
and it is felt that no further explanation is required [Ref. 
5] . 


The process of application of the boundary conditions and 
transformation of the system to a form suitable for solution 
is performed by subroutine BCOND. In the case where all pre¬ 
defined displacements are equal to zero the process can be 
considerably simplified. In this case the final form of the 
system is: 
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It is apparent that the only operation required is the substi¬ 
tution of the diagonal element by B. 

For this reason subroutine BCOND was provided with two 
entries, one for the case of all zero predefined displace¬ 
ments and one (BC0ND2) for the case of one or more nonzero 

20 

predefined displacements. A number B=10 is used and is 
considered sufficient for the purpose of the program. The 
operations are carried out in a step by step procedure in a 
sense that at maximum two stiffness matrix coefficient blocks 
and their corresponding load vectors are simultaneously present 
in core. 
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6. Solution of the System - Petemination of Reactions 


The solution of the system is performed by subroutine 
SOLVE. The technique employed was devised by Professor Gilles 
Cantin of the Naval Postgraduate School and presented in Ref. 4. 

In the process of solving the system of equations the 
inversion of diagonal blocks of the stiffness matrix is re¬ 
quired. This inversion is performed by subroutine SYMINV. 

Before inversion of the block is performed it is checked 
whether it is singular or nearly singular. The trace of the 

matrix is evaluated (excluding diagonal elements containing the 
20 

number B=10 ) and all diagonal elements are compared to the 

-12 

trace multiplied by a factor of 10 . If a diagonal element 

is found smaller than this product the block is characterized 
as nearly singular. The execution continues but a proper indi¬ 
cation is given by SYMINV to SOLVE which in turn gives a warning 

message at the output via subroutine ERROR. The choice of the 

-12 . • • 
factor 10 was cased on,the number of significant digits 

provided by an IBM/360 machine on which the program was tested. 
Should the program be used on a different machine, alteration 
of this number should be considered in a sense that the dif¬ 
ference between the exponent of the factor and the humber of 
significant digits provided by the machine should be adequate 
to provide for loss of accuracy due to round off errors. An 
extensive presentation of the concepts underlying this idea 
is beyond the scope of this work. 
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If one diagonal element is found to be equal to 
zero the block is singular and further continuation of the 
execution is impossible. Again a proper indication is given 
by SYMINV to SOLVE which in turn gives an error message at 
the output via subroutine ERROR and the execution is stopped 
for this problem. Consequently the calling procedure is re¬ 
started from the beginning and the execution of the next 
problem begins (if such exists). 

Upon inversion of each diagonal block [a] of the stiff¬ 
ness matrix a condition number c is evaluated for this block 
within SYMINV. This condition number is defined as 

= I|a|| ♦ Ia whereby |1 a|| is meant a noiin of [A] . This condi¬ 
tion number is printed out, and provides an upper bound on the 
expected relative error of the results of the problem. The 
difference between the number of significant digits provided 
by the machine used and the maximum exponent of the condition 
number gives approximately the number of significant digits 
in the displacement results. Further information on the 
significance of the condition number can be found in Ref. 6. 

Following the solution of the system of equations the 
reactions are evaluated by subroutine DISP as explained in 
Section B5 but with an opposite sign. Consequently the nodal 
displacements and reactions are printed out by DISP. 

7. Evaluation of Nodal Strains and Stresses 

The nodal displacements as found by the solution of 
the system 2.5 are’ stored in vector form and in a sequence 
following the numbering convention of the nodes for the whole 
structure (Section II.A). 
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The determination of element nodal displacements by 
use of relation 2>2 requires a displacement vector containing 
the displacement components of all nodes of the corresponding 
element in a sequence following the numbering convention for 
individual elements (Section II.A). The formation of such 
vectors is accomplished by properly sorting the displacement 
components via subroutine DISK (entry WDISKB) called by DISP. 

Subroutine STRESS evaluates the element nodal strains 
and stresses by use of relations 2.2 and 2.3. The evaluation 
of the product [D]{u} requires the use of subroutine FORMK. 

The computational details of this operation are explained in 
Re f. 1. 

Printing of element strains and stresses is optional 
since they provide useful information only for nodes on inter¬ 
faces between elements of different materials. 

Evaluation and printing of average nodal strains and 
stresses is performed by subroutine MAINE and is the conclud¬ 
ing operation in the solution of each problem. 

C. SUBROUTINE DISK 

Transfer of data to and from the data files stored on 
direct access storage devices is accomplished by subroutine 
DISK. This siibroutine uses multiple entries, each entry per¬ 
forming a different operation. Calling of DISK has been reduced 
to a minimum to avoid useless expense of execution time in data 
transfer. The most critical entries are those used for trans¬ 
ferring element stiffness matrices and total stiffness matrix 
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blocks. The reason being the size of the arrays transferred 
and the number of transfers required. Reduction of the 
transfer time requires the use of the maximum possible record 
length. The size of element stiffness matrices is known for 
each type of element, thus the optimum record length is deter¬ 
mined and the corresponding entries are coded on this basis. 
However, the block size in the equation solver is determined 
by the user and thus the optimum record length must be deter¬ 
mined by the user (see Section IV). The function of the 
various entries is explained below. Entry names starting with 
the letter W denote transfer of data to a data file and names 
starting with the letter R denote transfer of data from data 
files to core. 


RDDISK, WRDISK: 


File number 7. Transfer of total stiffness 
matrix blocks. 


RDISKl, WDISKl; File number 8. Transfer of element stiff¬ 
ness matrices. 

RDISKC, WDISKC: File number 10. Transfer of nodal coordinate 

component vectors COORD (3) 


RDSKRl, WDSKRl: File number 20. Transfer of load vector 

blocks when stored in work vector RBI. 
Similarly for RDISKR2, RDISKR3, WDISKR2, 
WDISKR3. 


RDISKB, WDISKB: File number 21. Transfer of individual 

node displacement components vectors. 

The random access data file thus formed is 
used for the formation of element nodal 
displacement vectors. 

RDISKS, WDISKS: File number 12. Transfer of vector SSSSl 

(See Appendix A.) 
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IV. PROBLEM AND INPUT DATA PREPARATION 


In this section the various steps in modeling a problem 
to be solved and the preparation of input data for the com¬ 
puter program are described. 

A. STRUCTURE MODELING 

The first step in applying a finite element solution to a 
problem is the selection of an appropriate mesh. In many 
cases an extrapolation of the results will be required and 
hence more than one meshes will have to be selected. 

The mesh size selection does not follow any predetermined 
rules and will, in general, depend on the nature of the problem 
and the judgement of the analyst. 

An arbitrary numbering convention must be adopted for the 
nodal points of the entire structure (in distinction with 
the numbering convention for the element nodes shown in Section 
II). The relation of node numbers in the two numbering systems 
establishes the connectivity matrix described in Section IV.B.l. 

Though no restriction is imposed on the selection of the 
numbering convention for the entire structure, it has been 
shown in practice that an intelligent selection may consider¬ 
ably reduce the bandwidth of the total stiffness matrix. 
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B. INPUT DATA PREPARATION 


Another factor to be determined by the analyst is the 
size of the blocks of coefficients in the equation solver. 

The block size affects the execution time required for solu¬ 
tion of the system of equations by subroutine SOLVE. Some 
experimentation has been conducted on this subject by Cantin 
and is presented in Ref. 2. On the basis of this experimenta¬ 
tion it seems that the larger the block size the lesser the 
execution time is. However, the subject has not, as yet, 
been thoroughly investigated and no firm rules can be stated 
at present. TRISOP allows the use of any block size. How¬ 
ever, if a block size smaller than the size of the element 
stiffness matrix is selected, the dimensions of the various 
arrays of common block B3 may have to be changed. 

After selection of the block size the number of records 
per block NREC7 and record length of file number 7 can be 
determined. Maximum transfer efficiency requires that NREC7 
is minimum and thus the record length the maximum possible for 
the machine used. Subsequently the remaining input data can 
be determined. In the following a description of the non 
self-explanatory data is given. 

1. Connectivity Matrix 

The size of this matrix is NELxNPEL. Each row cor¬ 
responds to one element and represents the correspondence of 
node numbers in the two numbering systems, namely, the ele¬ 
ment node numbering system and the structure node numbering 
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system. The values of the elements of each row are the node 
nximbers in the entire structure system and are sequenced 
according to the element node numbering system. Thus, the 
second element of the fifth row represents node number two 
of finite element number five and its value is the number of 
this node in the structure system. 

2. Nodal Point Coordinates 

The nodal point coordinates are defined in the global 
system of reference. Cartesian or cylindrical coordinates 
may be used. These coordinate systems are defined in Fig. 

3. All three coordinates of each node must be defined in the 



X 


Figure 3. 
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same system. However, different nodes may be defined in 
different coordinate systems. Each set of nodal coordinates 
must be accompanied by a code number (KND) indicating whether 
cartesian or cylindrical coordinates are used. 

3. Material Properties 

The maximum permissible number of different materials 
is ten. In case the structure is composed of more than one 
materials a code number (NC0N21) must be assigned to each 
material (NC0N21=1,2,..., n, where n is the number of materi¬ 
als present in the structure). Material properties are given 
in the form of matrix ELDAT. Each row represents a different 
material and its elements are in sequence. Young's modulus, 
Poisson's ratio and coefficient of thermal expansion. The 
present form of TRISOP does not provide for thermal stress 
analysis and thus no value is required to be given to the co¬ 
efficient of thermal expansion. 

C. DATA CARDS 

In order to preserve the flexibility of the input no read¬ 
ing formats are specified in the program except for titles, 
object time format arrays and general characteristics of the 
problem. Object time formats, stored in the array FMT, are 
used for reading input data. The same formats are used for 
echo-printing of the input data. In the description of input 
data cards the variable FMT denotes the object time format 
under which the corresponding data wll be read. One data 
card is required per format. Formats may be arbitrarily 
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specified by the user. However it is recommended that only 
one data card be used for the variables described in Sections 
IV.C.S.c; 6.c; 7.c; 8.c; 9.c and 10.c. The number of cards 
required per group of data is based on the above recommendation. 

Every group of data cards is preceded by a card containing 
alphanumeric information. Exception is made only for the con¬ 
nectivity matrix of cubic elements where two such cards are 
required. The contents of these cards are stored in the array 
TITLE and appear as identification headings at the output. 

Both TITLE and FMT are read under a (10A8) format. 

1. Characteristics of Data Files 

These data are given only once for each run and are 
read under a (7110) format. Six cards are required on which 
the values of record length and number of records for the 
various data files are given in the sequence: NR7, LRW7, NREC7, 
iMR8, LRB8, LRW3, NR9, LRB9, NRlO, LP310, LRWIO , NR12, LRB12, 
LRW12, NR13, LRB13, LRW13, NR14, LRB14, NR15, LRB15, NR16, 
LRB16, NR17, LRB17, NR19, LRB19, NR20, LRB20, LRW20, NR21, 
LRB21, LRW21, NR22, LRB22. 

These cards can be produced by use of the auxiliary 
program JCL as described in Appendix B. The required values 
of the above variables are given in the same Appendix. 

2. Initiation Card 

Format: 1A8 

Number of cards per problem: one 

The word START must be punched starting at column one. 
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3. Problem Identification 


Format: (10A8) 

Number of cards per problem: one 

Any alphameric information identifying the problem 
and punched in columns 2-60. 

4• General Characteristics of the Problem 
Format: (1015) 

Nur±)er of cards per problem: one 
Columns Variable Meaning 


1-5 

NEL 

Total number of elements 

6-10 

NDPT 

Total number of nodal points 

11-15 

NS 

Stiffness matrix bloclc size 

16-20 

NPEL 

Niimber of nodes per element 

21-25 

NMAT 

Number of different materials 

26-30 

NCLD 

Total number of nodes with 
concentrated loads 

31-35 

NPBC 

Total number of nodes with 
boundary conditions 

36-40 

NPBC2 

Number of nodes with non zero 
predefined displacements 

41-45 

NLINE 

Number of lines to be printed 
per page for strain and stress 
tables 

46-50 

KEL 

KEL=0: No printing of element 

nodal strains and stresses is 
required. 



KELt^O: Printing of above is 

is required. 

5 . 

Connectivity Matrix 



a. TITLE 

b. FMT 

c. Connectivity data. 
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Number of cards: NEL 

Variables: NCON= Array of connectivities for one 

element (one row of connectivity 
matrix) 

NC0N21= Material code number for the 

corresponding element (if required) 
The connectivity matrix is read row-wise one row at a 
time. In the case of structures composed of more than one 
material, the last element of NCON must be followed by NC0N21. 

6. Nodal Coordinates 


a. TITLE 

b. FMT 


c. Coordinate data 
Number of cards: 
Variables; IEL= 


COORD= 


KND= 


NDPT 

Node number in the entire struc¬ 
ture numbering system. 

Array of coordinates of node lEL 
in the global system. Coordinates 
are given in the sequence x, y, z 
or r, <}), z . 

Code number indicating cartesian 
or cylindrical coordinates. 

KND=0: Cartesian coordinates 

KND^O: Cylindrical coordinates. 


7. Material Properties 


a. TITLE 

b. FMT 
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c. Material property data 

Number of cards: NMAT (maximum ten) 

Variables: N= Material code number 

ELDAT= Array of material properties. 

8. Nodal Consistent Load Data 

a. TITLE 

b. FMT 

c. Load data 

Number of cards: NCL= Node number in the numbering 

system of the entire struc¬ 
ture 

CL0AD1= Array of concentrated load 
components in the global 
system. Load components 
are given in the sequence 

x,y,z. 

Load data cards must be in ascending order of NCL 

values. 

9. Boundary Condition Codes 

a. TITLE 

b. FMT 

c. Codes 

Number of cards: NPBC 

Variable: NBC= Array of boundary condition codes 

for each node. The first element 
is the node number and the remain¬ 
ing three represent the boundary 
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conditions applied in the x,y,z 
direction respectively. 

The codes employed are: 

0: Free node. 

1: Node constrained to zero 

displacement. 

2; A predefined displacement 
exists. 

The boundary condition code data cards must be 
sequenced in ascending order of node number. 

10. Predefined Displacements (if required) 

a. TITLE 

b. FMT 

c. Displacement data. 

Number of cards: NPBC2 
Variables: NDBC = Node number. 

BOON = Array of predefined displacement 
components in the global system, 
given in the x,y,z sequence. 
Predefined displacement data cards must be sequenced in 
ascending order of the values of NDBC. 

In case no predefined displacements exist in a problem 
(NPBC2=0) 3-11 cards described in this paragraph must be omitted. 

11. Subsequent Problems - Stop Trap for Last Problem 
In case more than one problem is to be executed in 

one run, the data information described in paragraphs IV.C,2 
to IV.C.IO must be provided for each problem. The initiation 
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card (START) of each problem must immediately follow the last 
card of the previous problem. 

The last problem must be followed by the following cards 
which serve as a stop trap. 

a. One card with the word "START" starting in 
column one. 

b. One card containing any alphanumeric information 
(or blank) . 

c. One card with any negative integer punched in 
columns 1-5. 

In order to provide for a better understanding of the 
above documentation sample inputs for two problems, one utiliz¬ 
ing quadratic and one cubic elements are presented in Appendix C. 

D. DIRECT ACCESS DATA FILES REQUIRE^ENTS 

The method of specification for the direct access file 
requirements will vary from one installation to the other. 

The particular features of the different installations will 
affect the form of subroutine DISK and the JCL cards. In 
Appendix B these requirements are described for the IBM 360/67 
operating system. The space requirements for each file as 
described in Appendix C are those to be used in establishing 
the form of subroutine DISK and JCL cards for any machine. 
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V. SOLUTION OF A SIMPLE PROBLEM 


When the initial form of TRISOP, for quadratic elements, 
was coded by Professor Cantin, it was tested on simple prob¬ 
lems, and results were obtained and compared with classical 
solutions. 

Testing of the present development of TRISOP was based on 
these problems and no further experimentation was conducted. 

In order to show the function of the program, as well as 
some observations affecting the use of it, the solution of 
a simple problem is presented in this section. 

The problem selected is that usually presented in classical 
texts of strength of materials as a simply supported beam of 
uniform cross section subjected to a concentrated load in the 
middle of the span. 

The model used for this problem for the quadratic element 
formulation is shown in Fig. 4. The concentrated load is sub¬ 
stituted by a line load. Nodes 1, 2 and 3 (here arbitrarily 
numbered) are constrained to zero displacement in all direc¬ 
tions and nodes 4, 5 and 6 are allowed to move in the z 
direction only. The model for cubic elements is the same with 
the difference that the constrained nodes are four at each . 
end instead of three. 

The consistent loads corresponding to the line load 
applied were found as required by the program. Such loads 
appear only at the nodes of the midsection of the beam and 
are shown in Fig. 6. 
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Five different meshes were employed in order to show the 
variation of results with mesh size. The meshes chosen were 

2, 4 , 6 , 8 and 10 elements in the z direction for both 
quadratic and cubic elements. 

Numerical integration was carried by use of 2, 3 and 4 
Gauss points (in different runs) for quadratic elements and 

3, 4 and 5 points for cubic elements. Listings of the sub¬ 
routines used are shown in Appendix I. 

Integration with 3 and 4 Gauss points for quadratic ele¬ 
ments and 4 and 5 points for cubic elements gave essentially 
the same results for this problem. For problems using 
elements with curved boundaries the results will vary with 
the number of Gauss points used. The solution with quadratic 
elements using 2 integration points was numerically unstable 
and no results were obtained. Numerical results for the 
maximum displacement at the midsection are given in Tables I 
and II. The same results are shown graphically in Fig. 7 for 
both quadratic and cubic elements. In the same figure, the 
execution time for subroutine FSTF and SOLVE, which consume 
most of the execution time for each problem, is shown. Figure 
8 shows graphically the displacements and execution time ob¬ 
tained from the solution using cubic elements. 

A. ALTERNATE FORMULATION 

In order to observe some further features of the finite 
element solution, the same problem was solved formulated in a 
different way. 
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By use of the existing syininetry in the z direction, only 
half of the beam was considered as shown in Fig. 5. All 
nodes lying in the plane of symmetry were constrained to a 
zero displacement in the x direction. The presence of addi¬ 
tional constraints enabled the use of two integration points. 
Numerical values for the maximum displacement in the midsection 
are shown in Table III. The same results as well as execution 
time are shown in graphical form in Fig. 9. 

B. DISCUSSION OF RESULTS 

Observation of the results presented here reveals that 

mesh refinement leads to improved results which eventually 
will converge to a certain value. Use of less accurate numeri¬ 
cal integration by reducing the number of integration points 
leads to improved results and considerable reduction of execu¬ 
tion time. Also the use of existing symmetry and in general 
different formulations of the problem may prove to be beneficial. 

Results obtained by use of cubic elements are considerably 
better than those of quadratic elements. However the execution 
time is appreciably higher for cubic elements and it is a 
matter of judgement whether cubic or quadratic elements should 
be used for each specific problem. 

These are observations essentially made on the results of a 
single basically simple problem. Thus no firm rules re¬ 

garding the use of TRISOP can be established and further ex¬ 
perimentation with different problems has to be conducted for 
this purpose. 
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TABLE I 


Quadratic Elements - Maximum Displacements in x 10 

ELEM. 2 4 6 8 10 

CUB 3 

or 5.761 6.847 6.992 7.026 7.038 

CUB 4 


TABLE II 

Cubic Elements - Maximum Displacements in x 10 

ELEM. 2 4 6 8 10 

CUB 5 

or 7.045 7.055 7.062 7.067 7.070 

CUB 4 

CUB 3 7.045 7.067 7.068 7.069 7.072 


TABLE III 

Quadratic Elements - Alternate Formulation 

-4 

Max. Disp. in X 10 

ELEM. 246 

CUB 4 

or 5.761 6.847 6.992 

CUB 3 

CUB 2 7.191 7.066 7.053 
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Figure 5 
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Quadratic elements 



Numbers represent proportions of total load applied on 
each nodec. 


Figure 6 
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APPENDIX A 


C0I4PUTER PROGRAM VARIABLES 

In this Appendix the most important variables appearing 
in the programs are listed together with a brief description 
of their function. The variables of each group are sequenced 
in alphabetic order. Arrays are shown with the dimensions 
required. 

A. DATA FILES CHARACTERISTICS 

NRX : Number of records of file number "X" 

LRBX : Length of records of file number "X" measured in 
bytes 

o 

LRWX: : Length of records of file number "X" measured in 
single words. 

NREC7 : Number of records per total stiffness matrix block 

NRX, LRBX AND LRWX are used only for checking the adequacy of 
disk storage space requested by the JCL cards. 

NREC7 specified by the user, is used in subroutine DISK 
and controls the transfer of stiffness matrix blocks. 


B. VARIABLES APPEARING IN 'COMMON' STATEMENTS 
AKl(IiS,NS) *, Work space used primarily for storage of 

stiffness matrix blocks when transferred in 
core. Alternatively is used for storing 
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I 



BCON (3) 


CLOAD(4) 


COORD (3) 


COREL(NPEL,3) 


equivalenced arrays and thus reduce the 
space requirements. Similarly for AK2 and 
AK3. 

: Array of non zero predefined displacement com¬ 

ponents for one node. It is used only in 
entry BC0ND2 of subroutine BCOND. The complete 
array of non zero predefined displacements is 
stored on disk, and only the components for one 
node at a time are present in core. 

; Concentrated load data for one node. The first 
three elements represent the components of the 
load and the fourth the node number. It is 
used only in subroutine FLOAD. The complete 
array of concentrated load data is stored on 
disk and only data for one node at a time are 
present in core. 

: Array of global coordinates for one node. The 

complete array of nodal coordinates is stored 
on disk in node sequence following the number¬ 
ing convention for the whole structure. 
Coordinates of only one node at a time are 
present in core. 

: Array containing the global coordinates for all 

nodes of one element. These arrays are formed 
by SUBROUTINE SORT one at a time and subse¬ 
quently stored on disk. Only one such array 
is present in core every time. 
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BLAST (616) 


ELDAT{10,3) 


KEL 


MM 


NBC(4) 


NBCH 


NBCL 


NCLD 


NCON (NPEL) 


Matrix of elastic constants of elements 
formed by subroutine ELPROP. 

Array containing information on properties of 
materials composing the structure. It is used 
in ELPROP for the formation of BLAST. 

Index number indicating whether printing of 
element nodal strains and stresses is re¬ 
quired or not. 

Number of stiffness matrix blocks per row 
(see Ref. 2) 

Boundary condition information for one node. 

The first element represents the node number 
and the remaining the boundary condition codes. 
The complete set of such arrays is stored on 
disk in node number sequence and only one is 
present in core at a time. 

Highest node number at which a boundary 
condition is applied. 

Lowest node number at which a boundary condi¬ 
tion is applied. 

Total number of nodes at which a concentrated 
load is applied. 

Array of connectivities for one element. This 
array shows the correspondence of element node 
numbers in the two numbering systems employed, 
namely, element numbering system and entire 
structure numbering system. The values of 
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the elements of the array represent the node 

numbers in the structure numbering system and 

are sequenced according to the element number¬ 
ing system. Thus the third element of the 

array represents node number three of the ele¬ 
ment and its value represents the number of 

this node in the structure numbering system. 

The complete set of these arrays is stored on 

disk and only one at a time is brought into 

core. 

NC0N21 

; Indicates the material of the corresponding ele¬ 
ment and is used in the case the structure is 

composed of more than one material. The values 

of this Variable must conform with the material 

number given in array ELDAT. 

ncount 

: Number of coefficients per stiffness matrix 

block. 

NDF 

: Number of special degrees of freedom. In the 

present form of the program its value is 

constant and equal to three. However, it is 

used throughout the origin as a variable in 

order to provide for future developments of 

the programs as for example degeneration to 

two dimensional elements. 

NDFP 

; Constant equal to: NDF+1 

NDPT 

: Total number of nodal points in the structure. 

NEL 

; Total number of finite elements used in the 

problem. 
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NEQ 

NLINE 


NMAT 


NN 


NPBC 


NPBC2 


NPEL 

NREC7 

NS 

NST 

NSTF 

SSSS1(7) 


UJNT (3) 


Total nuinber of equations of system (2.5). 

Number of lines to be printed per page in the 
output for displacement, reactions, strains and 
stresses. 

Number of materials of which the structure is 
composed. 

Number of stiffness matrix blocks per column (See 
Ref. 2) . 

Total number of nodes at which boundary conditions 
are applied. 

Total number of nodes at which non zero predefined 
displacements are applied. 

Number of nodes per element. 

Number of records per stiffness matrix block. 
Stiffness matrix block size. 

Size of element stiffness matrix. 

Number of coefficients per element stiffness matrix. 
Work vector used for transferring to and from data 
files the algebraic sums of stress or strain 
components contributed by elements sharing a common 
node. The seventh element of the array is used to 
count the number of occurrences of the node. Such 
arrays are formed in subroutine STRESS and further 
used in MAINE for the evaluation of average 
strains and stresses. 

Array of displacement components for one node. The 
complete set of such arrays is stored on disk and 
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displacements for only one node at a time 
are present in core. 


C. OTHER IMPORTANT VARIABLES 


In order to avoid redundant definitions, the most important 
variables of each subroutine that have not been presented in 
Sections A and B of this Appendix are explained below. 

1. Subroutine ELPROP 

E : Young's modulus 

PR : Poisson's ratio 

2. Subroutine CUB 

STK(NST,NST) ; Element stiffness matrix 

X,Y,Z : Gauss points coordinates 


3. Subroutine 
AJ(NDF,NDF) : 

AK(NST,NST) ; 

B(6,NST) : 

CORDG(NPEL,NDF): 

DNX(NDF,NPEL) 


FORMK 

Jacobian matrix 

Integrant of stiffness matrix integral 
formula 

Matrix [N*] = [D] [N] 

Element nodal coordinate array in the local 
system 

Array of derivatives of shape functions with 
respect to x,y,z. 


W1 (NDF,NPEL) 


Array of derivatives of shape functions with 
respect to C/n , . 


4. Subroutine BCOND 

: Stiffness matrix row number corresponding to 

the 'z' component of the boundary condition 
examined. 
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IBL,IBH : 

INDEX 

AKWl(NCOUNT) : 

AK1(NS,NS) 

AKW2(N COUNT) : 

AK2 (NS,NS) 

RBl(NS) : 

RB2(NS) : 

MMl 

IFR : 

5. Subroutine 


Lowest and highest row number in the block 
examined. 

Number of boundary conditions applied. 
Diagonal stiffness matrix block. 

Offdiagonal stiffness matrix block. 

Load vector corresponding to row 'I'. 

Load vector corresponding to blocks symmetric 
to the offdiagonal blocks of row 'I'. 
Specifies the extreme offdiagonal 
block to be examined in each row in order 
to avoid blank blocks. 

Specifies whether any boundary conditions 
have to be applied to the diagonal block of 
row 'I' or not. 

SOLVE 


KMM ; Same as MMl in MERGE. 

IFLG : Indicates whether block AK2 is singular or 

nearly singular. 


6. Subroutine DISP 


JLL 


Indicates the number of non blank elements in 
the last load vector 


7. Subroutine STRESS 


CORDG (NPEL, 3) 


SNELM(NPEL,6) 

SSELM(NPEL,6) 


Vector of element nodal coordinates in the 
local system. 

Array of element nodal strains. 

Array of element nodal stresses. 
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SSNN (6) 


Strain vector for one node. 


SSSS(6) 
UEL(NST) 


Stress vector for one node. 

Vector of element nodal displacements. 


i 
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APPENDIX B 


DIRECT ACCESS DATA FILES REQUIREMENTS 


In this Appendix the direct access data file require¬ 
ments as well as the necessary control cards for use of the 
compiler of release 18 of an IBM 360 operating system are 
described. 

A. DATA FILE SPACE REQUIREMENTS 

The data file space requirements are described by the 
variables NRX, LRBX, LRWX as they are explained in Section A 
of Appendix A. 

The required values of these variables are as follows: 
(the values of LRWX are not presented since always LRWX= 
LRBX/4) 

NR7 = MM X NN X NREC7 

LRB7 = NCOUNT X 8/NREC7 

NR8 = NEL X 4 for quadratic elements 


NEL X 12 for cubic elements 


LRB8 


7200 for quadratic elements 


6144 for cubic elements 


NR9 


NEL 


LRB9 


80 for quadratic elements 


128 for cubic elements 


NRIO 


NDPT 


LRBIO 


2 4 
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NR12 

= 

NDPT 

LRB12 

= 

48 

NR13 

= 

NDPT 

LRB13 

= 

48 

NR14 

= 

NCLD 

LRB14 

= 

32 

NR15 

= 

NPBC 

LRB15 

= 

16 

NR16 

= 

NEL 

LRB16 

= 

4 

NR17 

= 

NPBC 

LRB17 

= 

32 

NR19 

= 

NEL 

LRB19 

= 

480 

NR20 

= 

NN 

LRB20 

= 

NS X 8 

NR21 


NN X NS/NDF 

LRB21 

= 

24 

NR22 

= 

NPBC2 

LRB22 

— 

24 


B. CONTROL CARDS 

Job control language cards for an IBM 360 operating 
system are of the form: 

GO.FTXXFOOl DD UNIT = SYSDA, SPACE = (LRBX, (NRX, 1) ) , DISP= 
(NEW,DELETE) 

where X and XX represent the number of the corresponding file- 
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These cards may be produced by the auxiliary computer program 
named-JCL and whose listing is presented in Appendix H. This 
program is capable of printing and punching the JCL cards as 
well as the data files characteristics (Section IV.C.l) for a 
series of problems in one run. One data card is required for 
each problem. This card must contain the variables: NEL, 

NPEL, MM, NS, NDPT, NPBC, NPBC2, NCLD, NREC7, NMAT punched 
in format (IX,1015). Following the data card of the last 
problem, a card with a negative integer punched in columns one 
to six right justified, must be given. The control cards re¬ 
quired for this program in order to be used in an IBM 360 OS 
are shown in Appendix F. 

In case the value of MM is not known, it can be computed 
by the program from connectivity data. In this case MM is 
given the value zero in the data card described above and 
which must be immediately followed by the data cards: 

1. Object Time Format (20A4 ) 

Number of cards: one 

The format under which connectivity data will be read 
is punched. This format must provide for reading NPEL integers 
at a time, representing one row of the connectivity matrix. 

2. Connectivity Data 

Number of cards: NEL x (number of cards per row of 
connectivity matrix as per specified format.) 

Each individual row of the connectivity matrix is 
punched in sequence under the specified format. 
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The input data of program JCL are printed in the out¬ 
put for checking purposes. A sample output is presented in 
Appendix E. 

C. DEFINE FILE STATEMENTS 

The DEFINE FILE statements contained in subroutine DISK 
and defining the space requirements for random access files 
must be adjusted to the needs of the larger problem of each 
run. The form of these statements is: 

DEFINE FILE X(NRX,LRWX,U,IX) 

where X represents the number of the corresponding file. 

These statements are also produced by program JCL as described 
in Section B of this Appendix. 
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APPENDIX C 


SAMPLE INPUT DATA DECKS FOR TRISQP 


A. QUADRATIC ELEMENTS 
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APPENDIX D 


SAMPLE OUTPUT FOR TRISOP 

A sample output for the solution of a problem by use 
of quadratic elements is presented in this Appendix. The 
problem chosen is that of a simply supported beam with a 
line load acting on the mid-section and discretized in a 
1x1x2 mesh. This problem is one of those discussed in Sec¬ 
tion V. 

The form of the output for a solution using cubic ele¬ 
ments will be identical except for the connectivity matrix. 
The form of the connectivity matrix for the same problem 
solved by use of cubic elements will be as follows: 
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SIMPLY S'JPPORTED BEAM 


cuaooatic element 


MESH 1X1X2 


NEL,. •.total number CF ELEMENTS. 2 

NOPT.,.TOTAL NUMBER OF NODAL POINTS. 32 

NMAT, • .NUMBER OF DTFFEt^ENT MATERIAL*^. 1 

(MAXIMUM IS 10) 

nS .PIOCK SIZE f^OP THE LARGE CAPACITY SOLVER.60 

NPBC. ..NUMEER OF NODAL POINTS WITH ROUND. COND. 6 

NP^'C-’. .numrer of nodal points with PRFDEF. DISPL. 0 

NCLD...NUMBER CF NODAL POINTS WITH CONC. LOAD. 3 


A B C 0 


F G H 1 J K L N Q 


P 


Q R 


S T 


CONN. MTR. 


6 4 ! 

IP 16 n 


2 

14 


3 5 8 7 11 Q 10 12 18 16 13 14 15 
16 17 20 19 23 21 22 24 30 28 25 26 27 


17 2C IS 
29 ’2 31 


NEU. ... TOTAL NUMb*^R OF EQUATIONS. ^6 

N^^AND. .HALF-RAND WIDTH OF THE SYSTEM. 60 

NN.NUMBER OF BLOCKS PER COLUMN. 2 

MM.NUMBER CF PLOCKS PER ROW. 2 

NfOUNT.NUMBER OF COEFFICIENTS REF: BLOCK. 3600 

NSTF. . .NUMBER OF CC^FFICIENTS PEP ELEMENT. 3600 


COORD. JOIN‘D X 


Y 


2 


KND 
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0.0 
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24. ('OOOC 
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24.00000 

0 

6 

2.00000 
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2.000CC 
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0 
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0.0 

18.00000 
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0.0 

2,00000 

12.00000 
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1.00000 

12.00000 
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15 
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12.00000 
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1 .ocooo 

2.00000 

12.00000 
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17 

1.00000 

0. 0 

12.00000 
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18 

2 .ocooo 

2.00000 

12.00000 
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lo 

2.00000 

1.00000 

12.OOOOC 
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2 .ocooo 
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I 2.COOOO 

0 

21 
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0. 0 
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2.00000 

6.00000 

0 
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0.0 

6.00000 

0 

25 

0.0 

2,00000 

0.0 

0 

26 

0.0 

1.00000 

0.0 

n 

27 

0.0 

o.c 

Q.G 

0 

28 

1 .coooo 

2.00000 

0. 0 

0 

29 

1.00000 

0.0 

O.C 

0 

30 

2 .coooo 

2. {<0000 

0.0 

0 

31 

2.00000 

1 .coooo 

0.0 

0 

32 

2.00000 

0.0 

0.0 

0 


P.RATIO ALPHA 


FIFMENT data 



1 30000000. 0.30000 0.0 
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»:,** STRESS *»*»>" 

TIME = 32.T8 SCCOmS(CPU = 3.08 S''CO‘Jr'S) 









APPENDIX E 


SAMPLE OUTPUTS FOR AUXILIARY PROGRAM JCL 

Two different sample outputs for the auxiliary program 
JCL are presented in this Appendix, one for a problem utiliz¬ 
ing quadratic elements and one for a problem utilizing cubic 
elements. Both problems are from those presented in Section 
V for a simply supported beam. 


73 



t 




i 



input data 


.total NUMPtP Cf ELEMENTS. 4 

NDPT. . .TO'^&L NUMBER OF NODAI POINTS. 56 

NMAT. .* number CF OIFFEPENT MATERIALS. 1 

(MAXIMUM IS 10) 

NS.FLOCK SIZE FOR THE LARGE CAPACITY SOLVER. 60 

NPPC...NUMBER OF NODAL POINTS WITH BOUND. CONO. 6 

NPBC2 . .number of NODAL POINTS WITH PREDEF. DISPL. 0 

NCLD. . .NU'^^PER OF NODAL POINTS WITH CGNC. LOAD. ? 

NN.NUMBER OF BLOCKS PER COLUMN. S 

MM.NUMBER CF ^LOCKS PER ROW. 2 

NOEL...number OF NODES PER ELEMENT. 20 

NREC7. .NUMBER OF PECCPDS OF FILE NO 7. 4 


JCL CARDS 


//GO.FT07F001 
//GU.FT :)3F001 
//GO.FT09F001 
//bO.FTlO'^COl 
//GO.FT12F001 
//G0.‘=T13F00I 
//G0.FT20F031 
//GO.FT21'^001 
//GL..FT14FC01 
//G0.FT15F001 
//Gn.FTl7F00l 
//GO.FT19F001 


00 UNIT = SYS0A,$PA^I=C:^7G0, O);’^^!) I »Oi^P = (Nr:W,(JELETE) 
00 urjIT = SYSnA,SOACF = l 7?.)0, (OOU ,1 )) ,i)IS^=(NEv;, DELETE) 
OD 11NIT = SYSDA, SPACt: = ( 0J8D, (0 00^ ,1)),DI<^P=( NEW,DELETE) 
DO UN IT = SYSDA,SPACED(0024, (0056,1)) ,DISo=(NEW,0 ELETE ) 
DD UNIT = SYSOA,SPACr=(0356,(0056,1) ) ,DISR= LNFW,DELET^) 
no UN IT = SYS DA,SPACE=( 0056, (0056,1 )),DI$P=(NErt,DELETE) 
no liNIT=SYSDA,SPACF=(0480,(000^,1)),DISP=(NEW,OELETE) 
DD UNIT=SYSDA,SRACE=(00 24, (0060,1)) ,DISP = (NEW,DELETE ) 
00 UNI'^=SYS0A,SPArE=(0C'^2 , ( 0 03:^ , 1 ) ) , D I S P= ( NF W , 0 ELET E ) 
DD UN IT=SYSDA, SPACE=(0016, (0006,1 ) ) , D I S P= ( NEvV , D EL ET E ) 
DO UN I T=SYSDA, SPArr= ( 03^2 , ( 003 g , 1 ) ) ,ni SP = ( NFW,DFI.^TE ) 
DD UNIT = SYSDA,SPACC=(04 80,(0004,1) ) ,DISP=(NEW,DELETE ) 


DEFINE FILE STATEMENTS 


OFFINF 

FILE "( 

24, 

ie'J0,U,M7) , R ( 

16, IPJO.U,!3), 

DSKOOl53 

1 

10( 

56, 

6,IJ, I 13 ) , 


DSK0U16O 

2 

12( 

56, 

1 4 , U , r 1 2 ) , 1 3 ( 

14,IJ,I13), 

DSKOOl70 


20( 

% 

120,U,I20),21( 

60, 6,U,I21) 

DSK03180 



GENERAL 

FILE DATA 






NR X = number 

CE -RECORDS OF 

11 f X 






LPBX = LENGTH 

UF RFO 

□FD OF 

FILt. X IN 

BY 

TE S 




LSwX = -11- 

-1 

1- 

-IT- 

WGROS 




NR 7 = 24 

NPE 

16 

NR9 = 

4 

NklO = 

5o 

NR12 = 

56 

LPBT =7200 

LRB3 = 

7200 

LRB^ = 

BO 

LRB10= 

24 

LRB12= 

56 

LPW7 =1400 
NREr“^= 4 

LPW& = 

IROO 



LEW10= 

6 

LRW12 = 

1 4 

NP13 = 56 

NR 1 4 = 

3 

NR 15 = 

6 

NR 16 = 

0 

NR17 = 

6 

LRB1?= 56 
LPW13= 14 

LRP 14 = 

32 

LRB15= 

It 

LPP16= 

0 

LRB17= 

3 2 

NRio = 4 

N4 2 0 = 

3 

NK21 = 

60 

Nt^22 = 

0 



LPR19= 440 

LPP ^0= 

480 

LRB21= 

24 

LF 322 = 

0 




LRVs20 = 

1?0 

LRV'21 = 

6 





MOTE ; IF DAT 

A OF FIL^S N1 

16 OP 2^ 

ARP 

Z^RD ■^HFS 

r FILES 

, APE NOT 

REQUIRED 
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INPUT DATA 


NEL,,..TOTAL NUMBER OF ELEMENTS. 8 

NO^T. .'.TOTAL NUMBE^ pp NOHAL POINTS. 172 

NMAT. . rNUMBEO OF DIFFERENT MATERIALS. 1 

(MAXIMUM IS 10) 

NS.BLOCK SIZE FOR THE LARGE CAPACITY SOLVER. 96 

NPBC...NUMBER CF NODAL POINTS WITH BOUND. COND. 8 

NP5C2..NUMBEP OF NODAL POINTS WITH PREDEF. OISPL. 0 

NCLD. . .NUMBER OR NODAL POINTS WI^H CONC. LOAD. 4 

NN.number CF BLOCKS PER COLUMN. 6 

MM.NUMBER OF BLOCKS PER ROW. 2 

NPFL. .. NUMBER OF NODES PEP ELEMENT. 32 

NREC7. .MUMPER OF RECORDS OF FILE NO 7. 12 


JCL CAROS 


//G0.FT07F001 
//GO.FTGBFOOl 
//GO.FTC9F001 
//GD.FTlOFOOl 
//G0.FT12F001 
//GC.FT13F001 
//GO.FT20F001 
//GO.FTpiEoOl 
//G0.FT14F001 
//GC.^TISFOOI 
//G0.FT17F001 
//GC.FT19F001 


DO UNIT =SYSDA,SPACE=(6144,(0144,1) ),0 ISP=(NEW,0ELETE) 
DD UNIT = SYSDA,SPACE = (6144, (0096,1) ),DIS?=(NEW,DELETE ) 
DO UNIT=SYSOA,SPACF=(3128,(0038,1)),DISP=(NEW,DELETE) 
DD UNIT =SYSDA,SRACE=(0024, (0 172,1) ) ,DISP = (NEW,DELETE) 
DO UNIT=SYSDA,SPACF=(0056,(0172,1)),DISP=(NcW,DELETE) 
DD UMT = SYSnA, SPACE = ( OC56 , ( 0 172,1 ) ) , D I $ P = ( NEW , DEL ETE ) 
DO UN1T=SYSDA,SPACF=(0768,(0006,l)),DI5P=(NEw,DELETE) 
DD UNIT = SYSDA,SPACE = (0024,(0 192,1)) ,DISP=(NEW,DEL FTE) 
DO UfUT^SYSDA,SPACC=(0C32,(0004,1)),DISP=(NEW,DELETE) 
DO UNIT=SYSDA,SPACE=(0016,(0008,1)),DlSP=(NEw,D[LtTe) 
DD UNIT = SYSDA,SPArE-( 003 2, (0 038,1) ) ,01SP = (NEw , jELtTE) 
CD UNIT=SYSDA,SRAf p=(0483 ,(003=,1) ),0 I SP=(NEW,0ELSTF) 


define file STATE>^ENTS 


1 

2 

'X 


DEFINE FILE 7( 
10 ( 
12 ( 
20 ( 


144,15 36,U,I7) ,0( 96, 1536 ,U , I 8 ) i 

172, 6,U,I10), 

172, 14,U,I12),13( 172, 14,U,113) 

6, l‘^2,U,I20) ,21( 102, o,U,I21) 


DSKODloo 
DSK00200 
OSKC0210 
DSK0a220 


GENERAL FILE DATA «^^***:c 


MRX = NUMBER OF RECORDS OF FIL^ X 

LPPX = length OF RECORD OF FILE X IN BYTES 


LPWX = -11- 

-1 1- 

-II- 

NR7 = 144 
LRB7 =6144 
LFW7 =1536 
NREC7= 12 

NRR = 96 
LKB? =6144 
LPW8 =1536 

nr Q = 
LRB^ = 

NR13 = 172 
LFP13= 56 
LRW13= 14 

NP14 = 4 

LPP14= ^2 

NR15 = 
LP615= 

NR19 = 8 

LRblQ= 480 

NP23 = 6 
LRB20= 768 
LRW20= 1^2 

NR2 1 = 
LPF21= 
LPW21= 


NOTE : IF DATA OF FILES NO 16 OR 22 


WORDS 


p 

NPIO = 1 

72 

NifcU = 172 

12® 

Lr B10 = 

24 

L'’F.12= 5 b 


LRW1C= 

6 

LilWl? = 14 

8 

NR 16 = 

0 

NR17 = 8 

16 

L ^ B 1 6 = 

0 

LRcl7= 

192 

NC22 = 

0 


24 

6 

LPB22= 

0 


AR t 

ZERO THESE 

FILE 

S 4RE NOT P.EOUIP.EO 
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APPENDIX F 


DEC K STRUCTURE FOR COMPUTER PROGRAMS IN IBM 360 OS 
A. TRISOP - QUADRATIC ELEMENTS 
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cc •* 
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coo 

CO'-' 

II 

II UJ 
K-O 
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oa 

U-UO 


o 

UJ 

Q 

z 

< 

a: 

H- 

oi 

o 

U- 


r-ooco(\Jir. 

(i: CO o o o 

ocorooo 

_,oCoo 


z< 

on -It ■li' ^ ‘J' 



vOOjvO'^aioO'^^ 
inm—‘Ofococoajoo 
00000'4-'^C0 

ooooooooo 

OOOOOOOOO 


^1- -K ^ 'X* 


n II II II II 

UJUJUJ UJUJ 
ooo^w^o 
<<TC<<I 
CLa a a a 
LOOOOOOOOO 


11 n II II II u II II II 

U^- U' tU UJ U J UJ UJ u 1 lU 
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B. TRISOP - CUBIC ELEMENTS 
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C. AUXILIARY PROGRAil JCL 
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APPENDIX G 


COMPUTER LISTINGS FOR TRISOP 

In this Appendix two computer listings for TRISOP are 
presented, one for quadratic elements and one for cubic 
elements. 

The two listings are basically the same. The minor 
existing differences are: 

1. Core Space Requirements 

The core space requirements for cubic elements are 
higher than those for quadratic- Those requirements are re¬ 
flected in the dimensioning of arrays appearing in COMMON 
and DIMENSION statements. 

2 . Input of Connectivity Data 

The requirement for two headings and also reading of 
connectivity matrix in two parts for cubic elements dictates 
an alteration of the corresponding part of subroutine INPUT. 

3. Subroutine DISK 

Entries RDISKl and WDISKl differ in quadratic and 
cubic elements because optimization of transfer time for ele¬ 
ment stiffness matrices requires a different partitioning of 
these matrices into data file records for each case. Because 
of the similitude of the two programs the listing of TRISOP 
for cubic elements is not presented complete. Only the parts 
that differ from the listing for quadratic elements are shown. 
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Listing for subroutine CUB is not given in this Appendix 
since there is a choice of different subroutines each corres¬ 
ponding to integration by use of a different number of inte¬ 
gration points. Computer listings for these subroutines are 
presented in Appendix I. 
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•SAME AS FOR QUADRATIC ELEMENTS 
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SAME AS FOR QUADRATIC ELEMENTS 
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APPENDIX I 


LISTINGS OF INTEGRATION SUBROUTINES 
A, QUADRATIC SLEKSNTS - TWO GAUSS POINTS 
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D. CUBIC ELSI''ISNTS - THREE GAUSS POINTS 
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